Turning Point (3.1)

Background

The following is a demo of a turning point analysis (Jones, 2017) to determine the change from enrichment to dilution in 4 Iowa catchments. Each catchment has daily N concentration and water discharge data from 2015 through 2019 collected by Iowa State University.

Map developed by Dr. Peiyu Cao.

Import

The raw data is stored in an Excel file, with each catchment separated on a different sheet. This function imports the data for each site and calculates the area-normalized discharge (water yield). This is then stored in a list for future access.

# import data
site_fc <- function(n){
  
  # one catchment per sheet
  c1 <- read_excel("Data/Peiyu/Catchments_raw.xlsx", sheet = n)
  
  # rename columns
  c1 <- select(c1, 
               date = Date,
               Q_m3_day = `Inflow (m3/day)`,
               c_mg_L = `[NO3] in (mg/L)`,
               Interpolation) %>%
    filter(Interpolation == "F")  %>% # only keep raw data
    mutate(site = n,
           name = excel_sheets("Data/Peiyu/Catchments_raw.xlsx")[n],
           year = year(date),
           month = month(date),
           Q_norm_mm = case_when(
             name == "KS" ~ (Q_m3_day/3067200)*1000, # divide Q by drainage area in m2
             name == "RS" ~ (Q_m3_day/4806900)*1000, # and multiply by 1000
             name == "LP" ~ (Q_m3_day/2596500)*1000, # to convert to mm
             name == "WW" ~ (Q_m3_day/2188800)*1000
           )) %>%
    select(-Interpolation)
  
  return(c1)
  
}

site_list <- lapply(1:4, site_fc)

Turning Point

Once the data is imported, we made a custom function that separates the data by year. It then calculates the Parametric Quantile Regression (PQR) using the Ricker function as described in Jones, 2017.

The output of the PQR model returns the coefficients (a, b, c) for the equation. These coefficients can then be used to calculate the turning point for each year. The custom function is then applied to the 4 catchments.

Note: discharge is filtered to retain values above zero (in order to perform log functions). Higher discharge values are not yet filtered.

PQR_fc <- function(n){
  
  year <- c(2015:2019)
  name <- excel_sheets("Data/Peiyu/Catchments_raw.xlsx")[n]
  
  # make fc to filter data by year
  year_fc <- function(y){
    
    sub <- as.data.frame(site_list[[n]]) %>% 
      filter(year == y, Q_norm_mm > 0) %>% # log won't work on negative values
      select(date,
             x = Q_norm_mm,
             y = c_mg_L)
    
    return(sub)
    
  }
  
  # separate all years
  year_list <- lapply(year, year_fc)
  
  # make a new fc to apply PQR model for each year
  model_fc <- function(i){
    
    df <- as.data.frame(year_list[[i]])
    
    # plug equation into model 
    model_rq <- rq(y ~ log(x) + x, data = df)
    
    coef <- model_rq$coefficients
    
    return(coef)
    
  }
  
  # apply to all years, individually
  model_list <- lapply(1:length(year), model_fc)
  
  # calculate turning point based on model coefficients
  tp_fc <- function(i){
    
    a <- model_list[[i]][[1]]
    b <- model_list[[i]][[2]]
    c <- model_list[[i]][[3]]
    
    # Q in mm, area-normalized
    tp <- -b/c
    
    return(tp)
    
  }
  
  tp_list <- lapply(1:length(year), tp_fc)
  
  # the output will be one table with the tp values (by year)
  
  tp_vector <- c(tp_list[[1]],tp_list[[2]],tp_list[[3]],tp_list[[4]],tp_list[[5]])
  
  df <- as.data.frame(cbind(year, name))
  df <- df %>% cbind(tp_vector)
  colnames(df) <- c("year", "name", "tp_Q_mm_norm")
  
  return(df)
  
}

# all turning points by catchment-year data point
PQR_list <- lapply(1:4, PQR_fc)
PQR <- rbindlist(PQR_list)

Results

Turning points (tp) as area-normalized discharge:

Plots

When the custom PQR function is run with the data for the 4 catchments, some of the turning points have negative values. This can be further explored by visualizing the modeled curves.

Building off from the custom PQR function, a new one was made to automate the generation of the following plots.